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We derive a uniformly frustrated XY model that describes two-dimensional Josephson-junction 
arrays consisting of rotating Bose-Einstein condensates trapped by both a harmonic trap and a 
corotating deep optical lattice. The harmonic trap makes the coupling constant of the model have 
a nonuniform parabolic dependance. We study the ground state through Monte Carlo simulations 
in a wide range of the frustration parameter /, revealing a rich variety of vortex patterns. 
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Josephson-junction arrays (JJAs), a network of super- 
conducting islands, have attracted much interest because 
they are well-controlled systems to study nontrivial phase 
transitions as well as macroscopic quantum phase coher- 
ence [1]. The application of transverse magnetic fields 
to the superconducting JJA leads to realization of the 
uniformly frustrated XY model (UFXYM) 
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Here, J > denotes the coupling constant, 0j the phase 
of the superconducting node at a site j, and (jf) near 
neighbors. The bond variables Ajjt satisfy the constraint 
^2Ajj> — 27r/, where the summation is taken over the 
perimeter of a plaquette of the junctions and / is the 
magnetic flux (vortex) piecing the plaquette in units of 
the flux quantum. The vortices induces the frustration 
for the stable direction of the order parameter's phase 
at each site. The competition of two length scales — 
the mean separation of vortices and the period of un- 
derlying lattice — yields a rich variety of ground state 
structures, which depend on the rational or irrational 
number of / [2|, |3, ]% [5|, [y| . Also, the nature of the finite- 
temperature phase transition for nonzero / is still not 
fully elucidated, while for / = it is interpreted as the 
Berezinskii-Kosterlitz-Thouless (BKT) mechanism. For 
/ = 1/2, in particular, it remains controversial whether 
there are two distinct phase transitions associated with 
breaking of the continuous symmetry of U(l) gauge and 
discrete symmetry of Z2 chirality, closely connected with 
unbinding of kink-antikink pair excitation at Ising-type 
domain boundaries [7[. 

Cold atoms in a optical lattice (OL) provide an ideal 
testing ground for the study of many-body physics asso- 
ciated with the model Hamiltonian in condensed matter 
systems [8|. The advantage is that the microscopic pa- 
rameters of the periodic potential can be precisely con- 
trolled. The cold-atom analogs of JJAs have been real- 
ized in a one- dimensional (ID) OL [1, [IqI, where many 
Bose-Einstein condensates (BECs) are separated by po- 
tential barriers along the lattice direction. Also, it has 
suggested that BECs confined by a 2D OL can mimic the 
physics of 2D JJAs [Til] - Recently, thermally activated 



vortex formation, associated with the BKT mechanism, 
in such a 2D bosonic JJA was observed through the direct 
imaging of the density profile [l2[ . 



In this work, we investigate the rotation effect, anal- 
ogous to that of a magnetic field for superconductors, 
on the 2D JJAs consisting of an atomic BEC. A recent 
experiment by Tung et al demonstrated periodic pin- 
ning effects for vortices in a BEC by the rotating OL 
[Hj]. Several theories suggested rich phase diagrams of 
vortex states due to the interplay between the vortex- 
vortex interaction and the periodic pinning potential 
TBI . HR \l7 \. However, they considered them only 
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for a few values of the filling factor, the vortex num- 
ber per unit cell of the OL (frustration parameter /). 
Here, we consider BECs in a 2D deep OL, where the 
condensate fractions are well localized at the periodic 
potential minima to form a 2D JJA. The application of 
rotation to this system realizes the uniformly frustrated 
bosonic JJA [l8[ . The mapping into the UFXYM is help- 
ful to study the equilibrium vortex structure in a wide 
range of rotation frequency, because direct simulation of 
the Gross-Pitaevskii equation is a time-consuming work. 
Also, the model provide simple approach to explore finite- 
temperature effects, which could provides a new ground 
to verify unresolved problems in statistical physics de- 
scribed above. In this paper, we clarify the equilibrium 
vortex configuration in the rotating bosonic JJA using 
Monte Carlo simulations of the UFXYM in a wide range 
of the frustration parameter /. Since we treat explicitly 
the trapping potential in addition to the OL, the site- 
site couplings become nonuniform and a finite-size effect 
is expected. 



First, we derive the UFXYM to describe the rotat- 
ing bosonic JJA combined with the harmonic trap. The 
BECs in a deep 2D OL can be mapped onto the XY 
model, where the amplitude of the condensate wave func- 
tion is frozen at each site, but its phase is still a relevant 
variable Here, we make use of this formalism for the 
rotating system. The many-body Hamiltonian of bosons 
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in a rotating frame with frequency Q = Qz is 
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where m is the atomic mass and g = Anh a/m the cou- 
pling constant with s- wave scattering length a. The field 
operator if; obeys the bosonic commutation relations. 
Conservation of the total particle number is ensured by 
the chemical potential \i. The external potential consists 
of two parts V ex = Vho+VoL- a centrifugal- force- modified 
harmonic potential Vho — vn(uj\ — Vt 2 )r 2 /2 + muo 2 z 2 /2 
and a 2D OL Vol = Vo [sin 2 (irx/d) + sin 2 (jry / d)] with 
the square lattice geometry and the spatial periodicity 
d. The minima of the 2D OL are located at the points 
}d = (jxjy)d with integers j x and j y . 

We assume that the laser intensity is large enough to 
create many separated wells giving rise to a 2D array 
of condensates. Still, the small overlap between the wave 
functions of adjacent wells causes quantum tunneling and 
can be sufficient to ensure overall coherence of the sys- 
tem. If the energy due to interaction and rotation is small 
compared to the energy separation between the lowest 
and first excited band, the particles are confined to the 
lowest Wannier orbit als. Following the analogy of a Bloch 
electron in a magnetic field, we take the Wannier basis 

as ^(r) = djWj(r) exp (im/%) f* A(r') ■ dr' , where 
A = ft x r is the analog of the magnetic vector poten- 
tial, Wj (r) the Wannier wave function localized at the jth 
well, and dj the boson annihilation operator. The nor- 
malization condition J drw? (r)wy (r) = dyy implies the 

total number N = = J2j Ny 

With this basis, Eq. ([2]) leads to the Bose-Hubbard 
model in the rotating frame [19[ 
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where denotes a sum over nearest-neighbor 

sites and tyy = - J drw?(r) (-h 2 V 2 /2m + Vol) Wj'(r), 
Ej = J drw?(r) (-ft 2 V 2 /2ra + V ex - /j) Wj(r), and Uj = 
g J dr\wj(r)\ A represent the hopping matrix element, the 
energy offset of each lattice site, and the on-site en- 
ergy, respectively. The effect of rotation is described 
by Ayy = (m/h) f**' A(r / ) • dr' with the constraint 
J2 U c Ayy = 27r/, where the sum is taken around any 
unit cell of the 2D array. The constant / is the frus- 
tration parameter, being given by the average number of 
vortices per unit cell: / = 2Qd 2 / '/€, with quantum cir- 
culation n = h/m. The Hamiltonian ([3]) predicts novel 
vortex properties and fractal quantum Hall features of 
the strongly interacting lattice bosons [l{J [5o[. Other 



methods of creating this "effective" magnetic field have 
been discussed [21] . 

If the number of atoms per site is large (Nj ^> 1), the 
operator can be expressed in terms of its amplitude and 
phase, the amplitude being subsequently approximated 
by the c number as dj c± ^Nje 1,6 ^ . Then, Eq. (j3j) reduces 
to 

h = - J2 cos ft - °y + 4ij0 - E Vj ° 2 
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where we have used the phase representation Nj = Nj — 



id/dOy 6j = 6y and the notation Jyy = Uv ^/N } Ny. 
This reduction is valid when Jyy/N? < Uj fll|. 

The first term of Eq. (@]) corresponds to the UFXYM 
with spatially inhomogeneous nearest-neighbor coupling 
Jjj/. To neglect the other terms and to estimate Jjj/, 
the equilibrium form of Wj and Nj must be calculated. 
We assume that the equilibrium density is determined 
by minimizing the last c number term of Eq. (|4]), which 
is the dominant contribution of the ground-state en- 
ergy. Then, Ej + UjNj = and the third term may 
be neglected automatically. Next, we apply the ansatz 
Wj(r) — uo(x—j x d, y-jyd)vj(z) with the site- independent 
transverse part uq(x, y) = (^cr)~ 1 e~^ x +y ^ 2a and the 
site- dependent longitudinal part Vj{z) [22[. Since the 
atoms are tightly confined by 2D OL, the contribution 
arising from the two-body interactions is negligible for 
the estimation of uq(x, y) and the the variational param- 
eter a can be obtained easily. The longitudinal part is 
approximated by the inverted parabolic form Vj(z) 2 = 
(lij/g m Nj){l - z 2 /R 2 z -), with g 1B = #/2trt 2 , the local 
chemical potential /ij = m(u 2 — (] 2 )(j^ ax — j 2 — j 2 )d 2 /2, 
and the Thomas- Fermi radius R 2 ^ = 2/ij /muo 2 . Here, 
Nj = for |j| > j max because of the harmonic confine- 
ment. Using the normalization condition J Vj(z) 2 dz = 1 
and Nj = N, we can obtain 
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For a given Vo we evaluate the variational wave function 
uo(x,y) to obtain the optimized value of a. Through 
Eqs. ([5]) and (|6]) with this optimized a, the parameter 
values in Eq. (|4]) as well as Nj can be fixed. 

Under these formula we investigate the ground state 
of this system. Following the typical experimental con- 
ditions such as 87 Rb atoms used in JILA experiments 
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[H Hal, we use N = 6 x 10 5 and a = 5.29 nm. The 
frequencies of the trapping potential are set as uo± = 
11.5 x 2tt and uo z = 50 x 27r, which gives a_\_ =3.2 //m. 
The lattice spacing is set as d = 5 /im. 

We confirm that the obtained distribution Nj is quanti- 
tatively consistent with that obtained from the numerical 
solution of the 3D Gross-Pitaevskii equation; the particle 
number at the central well is iV( 0j o) — 6000, decreasing 
from the center to the outside according to Eq. (|6]). The 
conditions of the Josephson regime, Jjj/ /N? <C Uj and 
Jjj/ ^> /7j, are certainly satisfied. The former condi- 
tion is valid because of Nj ^> 1, even for outermost sites 
with Nj ~ 100. For the central region (j x ,jy) = (0,0) 
Uxify) ~ (1, 0) , the condition Jjj/ ^> Uj is well satisfied 
for 30 < V /huj_ < 90. We take V = 65huj ± in the 
following discussion, having J(o,o),(i,o)/^(o,o) — 100 and 
J(o,o),(i,o) = 0.9025ficj_L. Even for |j| ~ j max , the con- 
dition Jjy ^> Uj is still good. Therefore, the quantum 
correction arising from the third term of Eq. (|4|) may be 
neglected in our problem. 

We perform Monte Carlo simulations of the Hamilto- 
nian 
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The form of the coupling energy is 
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where we have used the optimized value of a and, 
when calculating the integral in Jjj/, the integral for 

R •/ 

the z direction was approximated as J_r., dzvjVj* ~ 

with Thomas-Fermi radius 
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R z j> and the area of the integral for the xy 
plane as dx dyuo(x, y)VoLUo(x, y). The symmet- 
ric gauge is chosen for the vector potential Aj ^. We 
use the Metropolis algorithm to study the ground-state 
properties of this system as a function of the frustration 
parameter /. For this purpose, the temperature is grad- 
ually decreased from high temperatures to zero accord- 
ing to the stimulated annealing. Since there are many 
metastable state caused by the frustration, we change 
the annealing rates in the several hundred simulations, 
taking the steady solution with the lowest energy as the 
ground state. 

It is known that the UFXYM of Eq. Q exhibits rich 
ground state structures depending on the parameter / [2|, 
y, 0]. For rational / = p/q, the ground state is periodic 
on the q x q cell in most cases. The striking difference of 
Eqs. (pp) and (|4j) of the bosonic JJA is the inhomogeneous 
coupling Jjj/ ex ^/NjNy. Also, it should be noted that 
the range of / is restricted by the harmonic potential 
because the rotation frequency £1 cannot exceed co± - 
that is, / < d 2 /7ra 2 L = 0.78 in our case. 




FIG. 1: (Color online) Ground-state energy and vortex lat- 
tice structures in a bosonic JJA under rotation. The top 
panel shows the total energy (normalized by J(o,o),(i,o) for 
/ = 0) as a function of /. The bottom panels from (a) to 
(f) represent the discretized condensate density iVy/ (black- 
white contour plot) and the positions of vortices marked by 
gray or red circles. Each square in the density corresponds 
to the site (minima of the OL), and vortices are located at 
the corners of the squares (maxima of the OL). The po- 
sitions of vortices are calculated by the current circulation 



9j/ + -Ay/) with the plaquette sum. The parameter 



values used are N = 6 x 10 5 , a = 5.29 nm, co± = 11.5 x 2tt 
Hz, u z — 50 x 2tt Hz, Vb = 65h(jj±, and a length of one side 
of the square, d = 5 /im. 



Figure [T] represents the total energy and the typical 
vortex patterns of the ground state as a function of /. 
The energy curve has a nonmonotonic behavior charac- 
terized by some minima at the simple rational values. 
These features are reflected in the bottom edge of Hofs- 
tadter butterfly spectrum [4]. The vortex configurations 
at these minima possess simple periodic structures as 
shown in Figs. []Ja)-(f), which represent the ground state 
for several values of / giving the visible minima of the 
energy curve. The vortices form a Bravais lattice with 
a unit cell of q x q and a quasi- ID structure oriented in 
parallel with one of the diagonals of the square lattice 
jjj 0, Uj]. This structure, called staircase states where 
constant currents flow along the diagonal staircases, was 
shown to be the true ground state for some limited values 
of / with simple rational forms such as / =1/2, 1/3, 2/5, 
3/7, 3/8 in the UFXYM with homogeneous coupling J 



4 



f= 2/9 - 0.222 f= 4/17= 0.235 /= 1/4 = 0.25 



/= 2/7 = 0.286 /= 3/10 = 0.3 /= 5/16 = 0.313 




FIG. 2: (Color online) The typical intermediate structures 
between / = 1/5 and / — 1/3. For / — 1/4 we also show 
two degenerate unit-cell structures of the ground state for the 
homogenous system. 



for the homogenous system as in Fig. [2] [2, these two 
configurations have exactly the same energy per site, and 
thus they are both ground states. Our simulations show 
that these two configurations are always separated by 
curved domain walls. In contrast, the variational result 
in Ref. [14( does not evidence the presence of degenerate 
configurations with the same energy. 

In conclusion, we derived a realistic UFXYM that de- 
scribes rotating BECs in both a trapping potential and 
a corotating deep OL. Monte Carlo simulations of this 
model clarify a variety of vortex phases for a wide range 
of the frustration parameter / that have not been pre- 
dicted by the Gross-Pitaevskii model. In future work, we 
plan to study finite-temperature properties such as an 
analog of competing phase transitions between the BKT 
type and the Ising type [3] in this inhomogeneous system. 

K.K. acknowledges supports of a Grant-in- Aid for Sci- 
entific Research from JSPS (Grant No. 18740213). 



While the periodicity of the vortex positions breaks 
slightly near the condensate edge, this staircase state can 
be the ground state for the inhomogeneous trapped sys- 
tem. For / = 1/2, a fully frustrated case, the vortex 
lattices form a checkerboard pattern, agre ement with the 
previous studies for trapped BECs PJ,[15L The energy is 
approximately reflection symmetric about / = 1/2 [23[, 
and the periodic structures for / > 1/2 are equivalent 
to those of 1 — /, but the condensate size is expanded 
and vortices are replaced by "vacancies"; an example is 
shown in Figs. [T] (d) and[T](f). 

Between these energy minima, we obtain characteris- 
tic intermediate structures consisting of the domains of 
simple periodic Bravais lattices; Fig. [2] shows an example 
of how one simple periodic structure (/ = 1/5) changes 
to another (/ = 1/3). Since the ground state has typi- 
cally q x q periodic unit cells, it is difficult to obtain the 
periodic structure for large q in the finite-size system. 
The periodicity is easily broken near the condensate edge 
due to the weak couplings [24] , the structural change be- 
ing of a crossover. This is contrast to the homogeneous 
model where the vortex patterns and accompanying do- 
main walls form diagonal lines for a square lattice, except 
for irrational values of / @ . This broken periodicity does 
not become noticeable as / increases, because the system 
size expands due to the centrifugal effect and approaches 
the homogeneous limit. For 1/3 < / < 1/2 the results re- 
produce the results obtained by the Coulomb gas model 
[5[. They consist of diagonal domains of the / = 1/2 
checkerboard configuration, separated by domain walls 
(or domains) of / = 1/3 structure. For / > 0.425, the 
ground-state structures are the / = 1/2 checkerboard 
pattern with a low concentration of missing vortices. 

An interesting case is for / = 1/4, where two possible 
vortex configurations of the ground state were proposed 
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